Help from code review

For the purposes of the code review I think we can ignore the mathematical part and focus on some other aspects:

• Check style guide ?
• Thoughts about the package as a whole?
• Documentation ?
• Any plots you think would be interesting
• Improvements on variable names ?
• Thoughts: the glmmGene is really legacy code since you can run glmmSeq with just one gene, should I remove it? Or just leave it out of the documentation ?
• Thoughts: if we allow further flexibility in the model then the plotting functions will not always apply. Should we do this and just add in catch errors is people try to use the functions incorrectly? Or is it better to keep the package simple for now and stick to our intended model use?

Outstanding tasks for developers

• add more commenting so its clearer to others what I am doing
• allow glm models (without random effects)
• allow model plot with a third variable (dummy)

Purpose

The aim of this package is to model gene expression with a general linear mixed model (glmm). This allows us to include random effects as well as mixed effects. For the purpose of the package we use the glmer function in lme4 package which fits glmm via REML or maximum likelihood.

Differences between other packages

Installing

To load the package as a library you can run the code in ./package_creator then:

library(glmmSeq)

Else you can source the functions individually:

functions = list.files(paste0(path_package, "/R"), full.names = TRUE)
invisible(lapply(functions, source))

But you will need to load in the additional libraries then:

library(MASS)
library(car)
library(ggplot2)
library(ggpubr)
library(lme4)
library(methods)
library(parallel)
library(plotly)
library(qvalue)
library(stats)
library(gghalves)

Overview

To get started first we load in the package

library(glmmSeq)

Then we can load in the example data set:

data(PEAC_minimal_load)
data(PEAC_dispersions)

metadata$EULAR_bin  = NA
metadata$EULAR_bin[metadata$EULAR_6m %in% c("Good responder", "Moderate responder" )] = "responder"
metadata$EULAR_bin[metadata$EULAR_6m %in% c("Non responder")] = "non_responder"
metadata = metadata[! is.na(metadata$EULAR_bin), ]

tpm = tpm[, metadata$SAMID]

The metadata contains the following, where EULAR is a rheumatoid arthritis response metric:

kable(head(metadata), row.names = F) %>% kable_styling()

Individual Genes

To model one gene over time using the model:

gene expression ~ fixed effects + random effects

In R the formula is defined by both the fixed-effects and random-effects part of the model, with the response on the left of a ~ operator and the terms, separated by + operators, on the right. Random-effects terms are distinguished by vertical bars (“|”) separating expressions for design matrices from grouping factors.

To achieve this we use a negative binomial model with gene dispersions passed in. These dispersions were previously created using DESeq2

Modelling a gene

For a specific gene, say MS4A1, we view the lmer fit using glmmGene:

metadata$Dummy = factor(sample(c("a", "b", "c", "d"), nrow(metadata), replace=T))

MS4A1fit <- glmmGene(~ Timepoint * EULAR_6m + Dummy + (1 | PATID),
                     gene = "MS4A1",
                     id = "PATID",
                     countdata = tpm,
                     metadata = metadata,
                     dispersion = disp["MS4A1"], 
                     pairedOnly=TRUE)

MS4A1fit

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.5589)  ( log )
## Formula: 
## count ~ Timepoint + EULAR_6m + Dummy + (1 | PATID) + Timepoint:EULAR_6m
##    Data: data
##  Offset: offset
##       AIC       BIC    logLik  deviance  df.resid 
##  619.0965  646.0944 -298.5483  597.0965        75 
## Random effects:
##  Groups Name        Std.Dev.
##  PATID  (Intercept) 0.6433  
## Number of obs: 86, groups:  PATID, 43
## Fixed Effects:
##                           (Intercept)                             Timepoint6  
##                               2.79223                                0.04196  
##            EULAR_6mModerate responder                  EULAR_6mNon responder  
##                               0.43192                               -0.97788  
##                                Dummyb                                 Dummyc  
##                              -1.04272                                0.23035  
##                                Dummyd  Timepoint6:EULAR_6mModerate responder  
##                              -0.48935                               -1.28024  
##      Timepoint6:EULAR_6mNon responder  
##                               0.05075  
## optimizer (bobyqa) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

Or to create a GlmmSeq object, we can use glmmSeq, which contains the following slots:

metadata$EULAR_bin = metadata$EULAR_6m

MS4A1glmm <- glmmSeq(~ Timepoint * EULAR_bin + (1 | PATID),
                     id = "PATID",
                     countdata = tpm["MS4A1", ],
                     metadata = metadata,
                     dispersion = disp["MS4A1"],
                     verbose=FALSE)

names(attributes(MS4A1glmm))

##  [1] "formula"        "stats"          "predict"        "reducedFormula"
##  [5] "countdata"      "metadata"       "modelData"      "optInfo"       
##  [9] "errors"         "variables"      "class"

Overall Model

To fit for all genes we can use the glmmSeq function. This creates an object of class GlmmSeq.

Note that this analysis can take some time:

• 100 genes takes a couple of seconds
• 1000 genes takes about 19 seconds
• All ~20000 genes takes about 9 mins

Therefore I will only run it for 1000 genes here.

genes = rownames(tpm)

resultTable <- glmmSeq(~ Timepoint * EULAR_bin + (1 | PATID),
                       id = "PATID",
                       countdata = tpm[genes, ],
                       metadata = metadata,
                       pairedOnly=TRUE,
                       dispersion = disp[genes])

## 
## n=86 samples, 43 individuals
## Time difference of 1.006739 mins
## Errors in 6 gene(s):IGHV1-58, IGHV1-67, IGHV3OR16-10, IL20RB, MMP7, THRIL

The variables used by the model are in the @modeldata slot

kable(resultTable@modelData) %>% kable_styling()

The model fit statistics can be viewed in the @stats slot:

stats = data.frame(resultTable@stats)
kable(stats[order(stats$P_Timepoint.EULAR_bin), ]) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

And the final model output in the @predict slot:

predict = data.frame(resultTable@predict)
kable(predict) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")